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Abstract: Joint modelling of longitudinal and survival data is inereasingly used in elinieal trials on 
eaneer. In prostate eaneer for example, these models permit to aeeount for the link between longitudinal 
measures of prostate-speeifie antigen (PSA) and the time of elinieal reeurrenee when studying the risk of 
relapse. In praetiee, multiple types of relapse may oeeur sueeessively. Distinguishing these transitions 
between health states would allow to evaluate, for example, how PSA trajeetory and elassieal eovariates 
impaet the risk of dying after a distant reeurrenee post-radiotherapy, or to prediet the risk of one spe- 
eifie type of elinieal reeurrenee post-radiotherapy, from the PSA history. In this eontext, we present a 
joint model for a longitudinal proeess and a multi-state proeess whieh is divided into two sub-models: a 
linear mixed sub-model for longitudinal data, and a multi-state sub-model with proportional hazards for 
transition times, both linked by shared random effeets. Parameters of this joint multi-state model are es¬ 
timated within the maximum likelihood framework using an EM algorithm eoupled to a quasi-Newton 
algorithm in ease of slow eonvergenee. It is implemented under R, by eombining and extending the 
mstate and 3M paekages. The estimation program is validated by simulations and applied on pooled 
data from two eohorts of men with loealized prostate eaneer and treated by radiotherapy. Thanks to 
the elassieal eovariates available at baseline and the PSA measurements eolleeted repeatedly during the 
follow-up, we are able to assess the biomarker’s trajeetory, define the risks of transitions between health 
states, and quantify the impaet of the PSA dynamies on eaeh transition intensity. 

Keywords: Joint modelling; Eongitudinal proeess; Multi-state proeess; Prostate eaneer; R; Shared ran¬ 
dom effeets. 
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1 Introduction 


In longitudinal health studies, marker data are usually collected at repeated measurement times until the 
occurrence of an event such as disease relapse or death, with the perspective to study the link between 
these two correlated processes or use the information brought by the marker’s dynamics to explain or 
predict the time to event. In such analyses, the repeated measurements of the marker should not be 


considered as a standard time-dependent covariate in a survival model (Andersen and Gill 1982 Fisher 


and Lin 19991 because the marker is an internal outcome measured with error and at discrete times 


whereas the Cox model assumes to know the exact values of the explanatory variables for all individuals 
at risk at each event time. To counteract these weaknesses, the two processes can be modelled jointly 
(Faucett and Thomas 1996t Wulfsohn and Tsiatis 19971. The principle is to define two sub-models 
(one mixed sub-model for the longitudinal data and one survival sub-model for the time-to-event data) 
and link them using a common latent structure. The shared random effect models, notably developed by 


Tsiatis and Davidian (2004i, are the most popular joint models. They usually assume that a function of 
the random effects from the linear mixed model is included as covariates in the survival model. Thus, 
longitudinal and survival processes share the same random effects, and the effect of the dynamics of the 
marker on the risk of event can be deduced. 

In prostate cancer, the joint modelling method is very useful. The prostate-specific antigen (PSA), 
which is a protein secreted by the prostate, is found to be over-expressed in the presence of prostate can¬ 
cer. This blood-based longitudinal tumour marker is commonly used by clinicians to monitor patients 
with localized prostate cancer following treatment (radiation therapy or surgery) in order to detect sub- 
clinical presence of disease. Proust-Lima et al. (20081, Taylor et al. ( 2005[ ) and Yu et al. (20081 showed, 
through various types of joint models, that the dynamics of this biomarker, along with the pre-treatment 
PSA level and other factors measuring the aggressiveness of cancer cells and the extent of the tumour, 
were risk factors for progression and permitted one to dynamically predict (i.e. using PSA to adapt 
prediction over time) the risk of clinical relapse. 


In practice, a patient may experience a succession of clinical progression events with for example 
a local recurrence, followed by a distant metastatic recurrence and then death. So, instead of the oc¬ 
currence of a single clinical event, the progression of prostate cancer should be defined as a multi-sfate 
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process with a focus on the transitions between clinical states and the impact of the biomarker dynamics 
on it. This is essential to understand and predict accurately the course of the disease, and it is of partic¬ 
ular relevance for the clinicians that need to distinguish the different types of events in order to properly 
adapt the treatment. 

Some authors already extended the classical framework of joint modelling to multiple time-to-event 
data. Chi and Ibrahim ( 2006|) proposed a joint model for multivariate longitudinal data and multivariate 


survival data. Liu and Huang (20091 and Kim et al. (2012 1 looked into the simultaneous study of three 


correlated outcomes: longitudinal data, times of recurrent events and time of terminal event. Elashoff 


et al. (20081 and Rizopoulos (2012|l extended the joint model to competing risks survival data, which 


allows to characterize the cause of survival event. Dantan et al. ( 2011| ) presented a joint model with 
latent state for longitudinal data and time-to-illness and time-to-death data. Recently, Andrinopoulou| 


et al. (20141 studied simultaneously two longitudinal markers and competing events. However, the joint 


study of longitudinal and multi-state data has never been proposed and implemented. Thus, we introduce 
a joint model with shared random effects for repeated measurements of a longitudinal marker and times 
of transitions between multiple states. It consists in a linear mixed model and a multi-state model with 
transition-specific proportional infensifies, bofh linked by shared random effecls. 

The compufafional aspecf is fhe main obsfacle in fhe developmenf of joinf models wifh shared ran¬ 


dom effecls. As explained by Gould el al. (2014|, fhe R package JM, developed by Rizopoulos (2010l, 
has enabled many advances in fhe use of join! modelling, particularly Ihrough efficienl numerical inle- 
gralions. On fhe olher hand, fhe R package instate, developed by |De Wreede ef ^ ( 2010[ l, provides 
esfimalion of mulfi-sfale models. In fhe presenf work, we combine and adapf fhese fwo packages in 
order lo esfimale join! mulfi-sfale models. Thus, fhe implemenlalion is easy and effeclive. Through fhe 
adapfalion of fhe j ointModel () funclion of fhe JM package, our approach uses fhe maximum likelihood 
approach, which is performed using fhe EM algorilhm coupled lo a quasi-Newfon algorilhm in case of 
slow convergence. The soflware advanlage is lhal if keeps fhe fealures, synlax and oulpuls of JM. 

The paper is organized as follows. Seclion presenls fhe joinf model for longifudinal and mulfi- 
sfale processes. The estimation and implemenlalion procedures are defailed in Seclion The joinf 
mulfi-sfale model has been validated by simulations, as defailed in Section]^ We apply our model on a 
defailed example of fwo cohorls of men wifh prosfale cancer in Seclion]^ A brief discussion is given in 
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Section |6] 


2 Joint multi-state model 


2.1 Notations 


For each individual i, a longitudinal process and a multi-state process are observed. Let t > 0) be 
the multi-state process where Ei(t) denotes the occupied state by subject i at time t and takes values in the 
finite state space S = {0,1,..., M}. It is assumed that the multi-state process is continuous and observed 
between the left truncation time r,o and the right censoring time C„ so that the observed process is 
Ei = {Ei(t), Tio < t < C,). We further consider that Ei is a non-homogeneous Markov process. The 
Markov property ensures that the future of the process depends only on the present state and not the 
past state, i.e. Pr(£',(t + u) = k\Ei{t) = h, {Ei{s), s < t}) = Pr(£',(f + u) = k\Ei(t) = h) ,Wh,k € S,Vm > 0 


(De Wreede et al. (2010 1 ), and the non-homogeneous property guarantees that the time since T/o impacts 
the future evolution of the process. Let us consider T, = (r,i, Ta ,..., the vector of the m,- > 1 

observed times for the individual i, with r,> < r,(r+i), Vr € {0,... ,m,' - 1), and where ^ denotes the 
transpose vector. If the last observed state for subject i (Ei{Ti„y)) is absorbing, that is it is impossible 
to leave it once entered (typically death), we observe m,- direct transitions. Otherwise, Ti^j equals Ci 
the right censoring time and we observe m,- - 1 direct transitions. We define by Si = (Sn ,..., the 
vector of observed transition indicators, with d,(,+i) equals 1 if a transition is observed at time Ti^r+i) 
(i.e. Ei{Tir) + £',(r,(r+i))) and 0 otherwise, Vr € {0,... ,m,- - 1). For each patient i, we also observe 
Yi = (Li,-the vector ofn,measuresofthemarkercollectedattimes ■■■, hn,, with < T™,- 


2.2 Joint multi-state model formulation 

The joint multi-state model is decomposed into two sub-models : a linear mixed sub-model for the 
longitudinal data (repeated measurements of the biomarker) and a multi-state model with proportional 
hazards for the event history data (transition and censoring times), both linked by shared random effects. 
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2.2.1 Longitudinal sub-model 

To model the trajeetory of the longitudinal marker, we use a linear mixed model. Under Gaussian as¬ 
sumptions, we assume that Yij the observed measure of the marker at time point is a noisy measure of 
the true measure Y*(tij). This non-observed measure Y*(tij) is explained aeeording to time and eovari- 
ates, at the population level, with fixed effeets yS and at the individual level, with random effeets bi that 
take into account the correlation between repeated measures of the same individual: 


Yij - Y*(tij) + eij 

- + + ( 1 ) 

with and Z,(t,y) the vectors of possibly time-dependent covariates associated respectively with 

the p-vector of fixed effecfs /3 and fhe ^-vector of random effecfs bi, bi ~ Af(0, D). Nofe fhaf e,- = 
{en,.. 6m, ~ Af(0, (T^Ini) where 1 is fhe idenfify mafrix; 6, and bi are independenf. 


2.2.2 Multi-state sub-model 

To model the transition times, we use a Markov multi-state model with proportional hazards that takes 
into account the marker’s dynamics through the shared random effects bi. Thus, for the transition from 
state h € S to state k € S, the transition intensity at time t takes the form: 


^ link 


Pr{Ei{t + dt) ^ k\Ei{t) ^ h-bi) 


dt^O dt 

s T, 


^ dhk.oit) ^^PiXhkj Jhk + Whk.iibi, f) r]hk). 


( 2 ) 


with Ahkfii-) the parametric baseline intensity (Weibull, piecewise constant or B-splines for example) and 
X^i^. the vector of prognostic factors associated with the r-vector of coefficients yuk- The multivariate 
function Whk.iibi, t) defines the dependence structure between the longitudinal and multi-state processes. 
We can choose Whk.iibi, t) = Y*it) (the true current level of the marker), or Whk.iibi, t) = dY*it)ldt (the 
true current slope), Whk,iibi,t) - (Y*it),dY*il)ldij (both), or any other function of the random effects 
in the context under study. Thus, the 5-vector of coefficients rjhk quantities the impact of the longitudinal 
marker’s dynamics on the transition intensity between states h and k. 
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3 Estimation 

3.1 Likelihood 

The parameters of this joint model are estimated in the maximum likelihood framework. Since the 
longitudinal and multi-state processes are independent conditionally on the random effects, the complete 
observed likelihood is obtained through the product of the individual contributions to the likelihood for 
the N individuals as: 



(3) 


where 6 is the vector of all the parameters contained in ([T]l and Q, and /(.) is a probability density 
function. 

In the longitudinal sub-part, described by the linear mixed model ([T]l, the conditional longitudinal 
outcomes are such that: 



(4) 


where ||x|| denotes the Euclidean norm of vector x, is the matrix of covariates with row vectors 
X^{tij)^, j - 1,...,«/, and likewise Z,- = {Z^hy)). 

For the multi-state side, let t) be the transition probability from state h to state k between the 
times s and t for individual i, i.e. = Pr(£',(0 - k\Ei{s) = h). For each r e {0, ...,m; - 1), 

the continuity and Markov assumptions imply that the individual i remains in state Ei{Tir) between the 
times Tir and Ti(r+\) with probability Ti(r+\)\bi), and if is an observed transition 

time, he transits to state £',(r,(r+i)) with intensity 4^ ^^, jj)(L(r+i)l^i)- By conditioning on Ei{TiQ), 

this translates in the individual contribution to the likelihood: 
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with 4^(0 = - HkMh possible delayed entry is here ignored by eonditioning on Ei(Tio). 

Linking the longitudinal and multi-state sub-parts, the random effeets bi follow a multivariate Gaus¬ 
sian distribution sueh that: 

3.2 Implementation 

The joint multi-state model has been implemented under R, via the eombination of two well-known 
paekages: mstate for the multi-state models and JM for the joint models with shared random effeets. 
To fit semi-parametrie Markov multi-state models, mstate eonsists first of preparing the database for 
multi-state analysis, more speeifically by defining eaeh pafienf’s hisfory as a series of rows, one for eaeh 
fransifion af risk for eaeh individual (in confrasl wifh only one dafa record (row) per individual in a 
classical survival analysis). By slralifying on fhe fransifion fype, fhe sfandard coxphC) funcfion of fhe 
R package survival can be used fo fif fransifion-specific Cox models. Wifh sfandard longifudinal and 
fime-fo-evenf dafa, fhe JM package inifialises fhe values of paramefers from fhe funcfion Imef) (nlme 
package) for fhe longifudinal sub-parl and coxph() (survival package) for fhe survival sub-parf. Then, 
fhe funcfion jointModelO carries ouf fhe estimation procedure. 

So by replacing fhe sfandard call fo coxph() by fhe call fo coxphf) on prepared dafa by mstate, an 
exfended jointModelO funcfion, called JMstateModelO, carries ouf fhe estimation procedure. The 
implemenfafion procedure is fhus disfinguished in four sfeps: 

• Ime O function fo inifialise fhe paramefers in fhe longifudinal sub-model; 

• mstate package fo adapf fhe dafa fo fhe mulfi-sfafe framework; 

• coxphO funcfion, on fhe prepared dafa, fo inifialise fhe parameters in fhe mulfi-sfafe sub-model; 

• JMstateModel () function fo estimate all fhe parameters of fhe join! mulli-slale model. 

A defailed example is given in Appendix C, and full defailed examples are available on https:// 
github.com/LoicFerrer/JMstateModel/. 



3.3 Algorithm 


The JMstateModelC) function computes and maximises the likelihood adapted to the multi-state data 
using integration and optimisation algorithms available in the JM package. Thus, the procedure combines 
an EM algorithm coupled to a quasi-Newton algorithm if the convergence is not achieved. Furthermore, 
the integral with respect to time in Q and the integral with respect to the random effects in Q do not 
have an analytical solution. These integrals are approached by numerical integration. The integrals 
over time are approximated using Gauss-Kronrod quadratures, and the integrals over the random effects 
using pseudo-adaptative Gauss-Hermite quadratures. Inference is provided by asymptotic properties for 
maximum likelihood estimators; the variance-covariance matrix of the parameter estimates is based on 
the inverse of the Hessian matrix. More details on the optimisation procedure, the EM algorithm and the 


numerical integrations are available in Rizopoulos (2012 1 . 


4 Simulation study 

The simulation study aimed to validate the estimation program described previously. 

4.1 Data generation 

For each subject i = 1,..., 1500 of the 500 replicates, the longitudinal and multi-state data were gener¬ 
ated according to the joint multi-state model defined as: 

‘ - iPo +/?o,x2f( + bio) + (/?i +/3i,x2fi + bn) x tij + eij, (7) 

'^‘hk^t\bi) = Ahkfiit) exp {yhkXi + t/m, level Y*(t) + r]hk,siopedY*(t)/dt), 

where the multi-state process that included three states (h, k € {0,1,2}) and three transitions is described 
in Figure [T] 


[Figure 1 about here.] 
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First, Xi and bi were generated according to normal distributions with respectively mean 2.04 and 


variance 0.5, and mean vector 


( \ 
0 

and variance-covariance matrix 

0.35 

-0.04 

0, 


-0.04 

V 

0.06 ^ 


. The times of mea¬ 


surements were generated such as Oy = 0,0.33,0.67,..., 16.33, and e,y was generated from a stan¬ 
dard normal distribution. The log baseline intensities were linear combinations of cubic B-splines 
with the same knot vector (0.004,4.120,7.455,10.908,18.201)^ for the three transitions, and the vec¬ 
tors of spline coefficients (-9.200, -3.500, -5.000, -3.900, -3.500, -2.500, -2.000)^ for the transition 
0 —> 1, (-9.860,-4.472,-5.128,-3.486, -2.457,-0.989,-0.715)^ for the transition 0 —> 2, and 
(-2.527, -2.170, -2.492, -2.156, -1.228, -0.955, -0.161)^ for the transition 1^2. Parameters values 
and knot locations were chosen according to the application data described in Section]^ 

The procedure to generate the vector of observed times T,- = (T/i,..., was inspired by 


Beyersmann et al. (20111 and Crowther and Lambert (20131. For each individual i, the censoring 
time C, was generated from an uniform distribution on [1,25], and the vector of true transition times 
T*q 2 , was generated according to the following procedure: (1) three random numbers 

M/.oi. Wi ,02 and were generated from three independent standard uniform distributions; (2) and 
r*02 were generated by solving “ Ag/^(voklbi) dvo/t + log(M,- o/t) = 0, for k = 1,2, through the Brent’s 


univariate root-finding method (Brent 19731; (3) then, the true transition time T*,, was generated via 


T* 

A\Jv\ 2 \bi) dvi 2 -i-log(M,- 12 ) = 0. Finally, by comparing T* and Q, the vector Ti, which characterizes 
the multi-state process, was deduced. 

The longitudinal measurements, generated from the linear mixed sub-model, were truncated at Tn 
the first observed time of the multi-state process. 


4.2 Estimated model 

The estimation model was the model defined in (|7|) wifh bi ~ N 


0 


vO, 

V / 


Dll 7^12 
Di 2 D22 




and eij ~ A((0, cr^). 


The log baseline infensifies were approximafed oy a linear combination of cubic-splines wifh fhree 
infernal knofs placed al Ihe quanliles of Ihe observed times. 
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4.3 Simulation results 

The simulations results were obtained through 500 replicates of 1500 individuals. Each joint multi¬ 
state model was estimated using 3, 9 and 15 pseudo-adaptative Gauss-Hermite quadrature points. The 
simulation results are presented in Table 

These results were very satisfying with unbiased estimates and correct 95% coverage rates. They 
showed however the need to use a certain number of Gauss-Hermite quadrature points to approximate the 
integral over the random effects. Indeed, the use of 3 pseudo-adaptative Gauss-Hermite quadrature points 
induced a bias in the estimation of the parameters associated with time effect in the longitudinal sub-part. 
This error was fully corrected with 9 points. Overall, these results confirmed the good performances of 
the implemented procedure. Complementary simulation results based on 500 and 1000 subjects are 
provided in Appendix A. 

5 Application 

We analysed data from patients with localized prostate cancer and treated by external beam radiotherapy. 
The analysis aimed to explore the link between PSA dynamics and transition intensities between clinical 
states, as well as to describe PSA repeated measurements and times of transitions between health states. 

5.1 Data description 

Our study focuses on 1474 men with a clinically localized prostate cancer and treated by external beam 
radiotherapy (EBRT): 629 patients come from the multi-center clinical trial RTOG 9406 (Radiation 
Therapy Oncology Group, USA) in which data collection has been conducted from 1994 to 2013, and 
845 patients come from the cohort of the British Columbia Cancer Agency (BCCA) in Vancouver, 
Canada, and have been examined between 1994 and 2012 (Table|^. During his follow-up, a patient can 
possibly go through several states defined as local recurrence, distanf recurrence, inifiafion of hormonal 
therapy (HT) and death, due or not to prostate cancer. The initiation of salvage hormonal therapy, 
which is an additional treatment prompted by physician observed signs in PSA or clinical signs, is 
designed to prevent growth of potentially present subclinical cancer. This intervention is not planned 
at diagnosis or initiated by any precise rule, but is rather based on a mutual agreement between the 
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clinician and his patient. Thus, it is treated as a disease state transition representing failure of the initial 


treatment to satisfactorily control the disease. Furthermore, as recommended in the article by Proust 


Lima et al. (20081, we only considered the local relapses which took place three years or later after 


radiation, or within three years of EBRT when the last PSA value was > 2 ng/ml. PSA data were 
collected at regular visits, for a median number of 10 PSA measurements per patient. Note that this 
number includes only PSA data recorded between the end of EBRT and before the first event (first 
clinical recurrence, hormonal therapy, death or censorship). Subjects with only one PSA measure were 
excluded, and subjects who had an event in the first year after EBRT were excluded to prevent inclusion 
of patients who had substantial residual initial tumors. As shown in Table three baseline factors 
were considered: the pre-therapy level of PSA in the log scale (iPSA), the T-stage category which 
characterizes the tumour size (3 categories were considered: 2, 3-4 versus 1 in reference), and the 
Gleason score category which measures the aggressiveness of cancer cells (3 categories: 7, 8-10 versus 
2-6 in reference). In the models, a cohort covariate was also considered coded as 1 for RTOG 9406 and 


-1 for BCCA. More details on the RTOG 9406 trial can be found in Michalski et al. (20051, and on the 


BCCA cohort in Pickles et al. (20031. 


[Table 1 about here.] 


The PSA individual trajectories collected between the end of EBRT and the occurrence of the first 
event are depicted in Eigure Overall, this longitudinal process is diphasic, with a decrease in the 
level of PSA in the first years following the end of EBRT, and a subsequent stabilisation or linear rise 
thereafter. According to the type of first relapse, the biomarker’s long-term increase may have different 
intensities (see “Hormonal Therapy” and “Censorship” for example). 

[Eigure 2 about here.] 

The multi-state data are depicted through the transitions between the 5 states and the corresponding 
amount of observed direct transitions in Eigure 

[Eigure 3 about here.] 

Erom the end of EBRT (state 0), a patient can experience either a transition to a localized recurrence 
(state 1), an hormonal therapy (state 2), a distant recurrence (state 3) or death (absorbing state 4). After 












12 


a localized recurrence (state 1), a patient may experience either a distant recurrence (state 3) or die (state 
4) or initiate a HT (state 2). After initiation of HT, a patient may only experience a distant recurrence 
or die, and finally, after a distant recurrence, a patient may only die. In total, 144 subjects had a local 
recurrence; 317 men initiated an hormonal therapy including 90 after a local recurrence; 90 men had 
a distant recurrence including 10 directly after a local recurrence and 33 after a HT initiation. In total, 
802 patients died including 523 who did not have another recorded progression of the cancer before. 
Among the 672 men who were censored during the follow-up, 533 were censored before experiencing 
any clinical progression. 

5.2 Specification of the joint model 

The joint multi-state model being a complex model, a step-by-step procedure was carried out to specify 
the joint model. The specifications of the longitudinal and multi-state sub-models were based on two 
separate analyses, that is assuming independence between the two processes. Covariate selection was 
made using uni- or multivariate Wald tests. 


5.2.1 Longitudinal sub-model specification 


The biphasic shape of log-PS A was described in a linear mixed model with two functions of time accord¬ 
ing to previous works (Proust-Lima et al. 20081: f\{t) = (1 -i- 0“ - 1 and / 2 (f) = /(I + tY, where 


a and v were estimated by profile likelihood (a = -1.2, v = 0). Thus, these two functions depicted 
respectively the short term drop in the level of log-PSA after EBRT and the long term linear increase 
of log-PSA. By denoting T,y = log ^PSA,(6y) - 1 - 0.1^ the log-measure of PSA for the individual i at time 
tij -the natural logarithm transformation is performed to obtain a Gaussian shape for the longitudinal 
response- the linear mixed sub-model took the form: 


= ipo + ^/3o,cov + ^io) + 

[Yi + ^Pl.cov + ^il) X flYij) + 

(/?2 + ^/32,cov + X flYij) + €ij, 
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with bi - {bio, bii,bi 2 y ~ N (0, D), D unstructured, and e,- = {en ,~ N{0, cr^Im)- The covari¬ 
ates Xp and Xp were sub-vectors of the baseline prognostic factors obtained using a backward 
stepwise procedure. For the sake of brevity, we will speak about PSA dynamics and biomarker’s cur¬ 
rent level/slope when referring actually to respectively the dynamics of log(PSA -i- 0.1) and the current 
level/slope of Y*{t). 

5 . 2.2 Multi-state sub-model specification 

In the multi-state sub-part, the determination of prognostic factors and proportionality between baseline 
intensities was also made by considering no link between the two processes (77 = 0) and unspecified base¬ 
line intensities (i.e. using a standard semi-parametric multi-state model). The full sub-model considered 
transition-specific baseline infensifies and fransilion-specific effecls of baseline prognostic faclors. To 
reduce fhe excessive number of paramefers fo esfimafe, a firsl sfep was fo assume proporfional baseline 
infensifies for some fransifions. Clinically, if made sense fo consider proporfional baseline intensifies for 
fransifionsleadingfolocalrecurrenceorhormonalfherapy: /loi,o(0 = exp(^02)To2,o(0 - exp(^i2)/li2,o(0; 
and for fhe fransifions leading fo disfanf recurrence: dos.oCO - exp(^i3)/li3_o(0 - exp(^23)T23,o(0- These 
assumpfions were confirmed by fhe dafa. We could nol make fhe same assumpfion for all fransifions 
leading fo deafh because fhe proporfional hazards assumpfion was nof verified. Instead, we chose 
Ti4,o( 0 - exp(^24)T24,o(0 and /lo 4 ,o(t) was slralified on fhe cohorf. This procedure reduced fhe num¬ 
ber of baseline intensifies fo six. A second sfep consisfed in selecfing fhe prognosfic facfors. Factors 
wifh an associated p-value > 0.5 were removed, and common covariafe effecfs on several fransifions 
were considered using mulfivariafe Wald fesfs. For example, fhe baseline T-sfage category had fhe same 
effecl on fransifion infensifies 0 ^ 1, 0 ^ 3 and 2^3. Finally, covariafes wifh p-value <0.1 were 
selecfed by using a backward stepwise procedure. 

5.2.3 Joint multi-state model specification 

In the joint model, log baseline intensities approximated by linear combinations of cubic B-splines with 
three internal knots replaced the unspecified ones. The dependence function Whk,i{bi, t) was the same 
for all the transitions h ^ k and was determined using Wald tests. It resulted that the combination of 
the true current level and the true current slope of the biomarker fitted at best the relationship between 



14 


PSA dynamics and the instantaneous risk to transit between health states. The final step was to select the 
prognostic factors and the log-coefficients of proportionality between baseline intensities with p-value < 
0.1 in the multi-state sub-part by using a backward stepwise procedure. Thus, the multi-state sub-model 
was: 

f f \\ 


Ak^Abi) = Tm.oCO exp 


^hk,i yhk + 


y*{t) 

dY*{t)ldt 


Vhk,level 

%/:,slope 


where the relations between Ahk,oit) and the final for /j, k e {0,..., 5) are indicated in Section 


5.2.2 


and in Table Note that the covariates that were removed of the joint model specification are not in 


Table [3l 


5.3 Results 

The parameter estimates of the joint multi-state model are presented in Table These parameters 
were those selected according to the procedure described previously. The parameters of the baseline 
intensities are not shown here for clarity. 


[Table 2 about here.] 

The estimated regression parameters in the longitudinal sub-part confirmed that pre-treatment PSA level 
was associated with the initial PSA level and the biphasic PSA trajectory, T-stage value was associated 
both with the short term and the long term dynamics while Gleason score was only associated with the 
long term trajectory. Higher values of these covariates measured at baseline corresponded to higher long 
term PSA levels. The cohort effect indicated a significant difference between the two cohorts only for 
the long term PSA evolution, with a greater long term increase of PSA in Vancouver. 

For the multi-state process, the prognostic factors showed that an advanced initial stage was not 
always associated with the intensities of transitions between health states after adjustment for the PSA 
dynamics. In particular, the Gleason score had significant effects on only two transition intensities. 
Moreover, we found that having a high PSA value at baseline was significantly associated with a higher 
instantaneous risk to directly experience hormonal therapy initiation or death after EBRT, but reduced the 
intensities of transitions leading to distant recurrence or death after a previous event. A poor (i.e. higher) 
T-stage category at baseline had globally a deleterious effect on the clinical endpoints. For the transitions 
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leading to distant recurrence after EBRT or after hormonal therapy, a patient with a Gleason score of 7 
at baseline had a 2.65 = exp(0.973) (95% Cl = 1.64^.28) higher hazard to transit than a patient with a 
Gleason score < 7. The cohort was significantly associated with the intensities of transitions leading to 
death after clinical recurrence or hormonal therapy -and the direct transition leading to local recurrence 
after EBRT. The instantaneous risk to experience these transitions was higher in BCCA. Eor the direct 
transitions leading to distant recurrence after local recurrence or hormonal therapy, the cohort effect was 
also significant, with higher intensities in RTOG 9406. 

Regarding the association parameters between PSA dynamics (current level and current slope) and 
clinical progressions, there were highly significant deleterious effects of the PSA dynamics from the 
initial state on the intensities of transitions leading to all the types of progression (local recurrence, hor¬ 
monal therapy or distant recurrence). Eor example, after adjustments for covariates and for the true slope 
of the biomarker, an increase of one unit of the true biomarker’s level (log PSA without error measure¬ 
ment) induced a 1.45 = exp(0.372) (95% Cl = 1.23-1.72) higher risk to experience a local recurrence. 
These results were expected: in patients with localized prostate cancer and treated by radiotherapy, a per¬ 
sistently high PSA level or/and a strong increase of PSA leads to higher hazard to experience a clinical 
recurrence or an additional therapy. In contrast, for the direct transition leading to death after radiother¬ 
apy, we found a deleterious effect of the current slope and a protective effect of the current level of the 
biomarker: at a given moment in the initial state, for two patients with the same baseline characteristics 
and the same slope of log PSA, the one with higher PSA value will be less likely to directly die. In 
this studied population, an important cause of death is induced by comorbidities, most of death from 
prostate cancer experienced a documented disease progression before. Prom the local recurrence, there 
was large deleterious effect of the current slope of the biomarker for the intensities of transitions leading 
to the hormonal therapy or the distant recurrence, and there was a borderline significant protective effect 
of the current level for the intensity of transition leading to the distant recurrence. Prom the hormonal 
therapy or the distant recurrence, there was no significant effect of the PSA dynamics on the hazard to 
change state. This was also clinically sensible, as it reflects that progression in these advanced stages is 
no more linked to PSA increase. In practice, criteria other than PSA are considered specifically in this 


phase of the disease, such as the PCWG2 criteria (Scher et al. 2008). Moreover, deaths in patients with 
hormonal therapy might be explained by cardiac toxicity due to HT. 
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5.4 Diagnostics 

The parameter estimates of the joint multi-state model were validated by several graphical tools pre¬ 
sented in Figure]^ For the longitudinal sub-model, the plotted standardized conditional residuals versus 
fitted values of the biomarker confirmed the homoscedasticity of the conditional errors. Subject-specific 
predictions were also compared to observations by plotting the average values by time intervals based 
on the deciles of the observation times. 95% confidence intervals of the observed values were added and 
confirmed the very good fit of the model concerning the longitudinal data. For the multi-state sub-part, 
we focused on P(0, t) - {PhkiO, t)) the matrix of transition probabilities between the times 0 and t. We 
compared our parametric estimator (obtained through the average of the predicted individual transition 
probabilities from the joint multi-state model) to the Aalen-Johansen estimator (non-parametric estima¬ 
tor of the transition probabilities), both using product integrals. This comparison is fully discussed and 
detailed in Appendix B. These comparisons showed the overall good performances of the joint multi¬ 
state model in terms of fit of transition probabilities, with the exception for the transition 1 —> 2 for 
which the immediate pike after EBRT could not be correctly captured by splines. 

[Figure 4 about here.] 


6 Discussion 

The joint model for the longitudinal biomarker PSA and multi-state clinical progression data provides 
a full model of prostate cancer progression which takes into account both classical prognostic factors 
and the dynamics of the PSA, in order to study factors that influence the transition intensities between 
clinical health states. The implementation is easy and relies on the mstate and JM packages. The multi¬ 
state data are prepared through the mstate package, and a slightly modified jointModelC) function 
carries out the estimation procedure. The estimation program has been validated by simulations, with 
very good performances. It underlined however some bias in the estimates when too few quadrature 
points (3) were used. In the absence of specific validation, we thus recommend at least 9 Gauss-Hermite 
quadrature points in the pseudo-adaptative numerical integration to approximate at best the integral over 
the random effects. Goodness-of-fit of the model can be assessed using diagnostic graphical tools whose 
methodology was detailed in Appendix B. 
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The application confirmed that the PSA dynamics strongly impacted the instantaneous risk to ex¬ 
perience a clinical recurrence or hormonal therapy after the end of radiotherapy. Moreover, the current 
slope of the biomarker had highly significant deleterious effect on the hazard to transit from local recur¬ 
rence to hormonal therapy or distant recurrence. Conversely, extrapolating the biomarker’s dynamics did 
not impact anymore the transition intensities from the hormonal therapy or the distant recurrence states. 
This expresses that in the advanced cancers, the PSA -and especially the collected measures prior the 
first event- is no more of importance. In these situations, other criteria have to be monitored. 

Previous works in prostate cancer had found a strong association between slope of log-PS A and any 
clinical recurrence (see Sene et al. ( |2014 i; Taylor et al. ( 2013| ), by considering all the recurrences in a 
composite event and the hormonal therapy as a time-dependent covariate. The limit of these approaches 
was that in practice the type of progression is of major importance as the care greatly depends on the 
type of risk the patient has. The joint multi-state model formalizes this need. In the same way as it 
was done with a single event (see Proust-Lima and Taylor (20091; [Rizopoulos (20111), individualized 
dynamic predictions of each type of progression could be derived from this model in order to precisely 
quantify the risk of each type of progression according to the PSA history. For example, the cumulative 
probability for subject i to reach the state k between the times s and t, s < t, given he was in state h 
at time s, could be expressed as: t) = £‘ Pr(Ei(u) = k\Ei{s) = h, xf ) du, with the 

history (i.e. collected measures) of the marker up to time s, X. the history of the longitudinal sub-part 
covariates until time s, and Xf = {X^^ .) the matrix containing the prognostic factors for all the state 
transitions. 

In this article, we assumed a continuous and Markov multi-state process. Depending on the topic, a 
semi-Markov process, which considers the spent time in the current state, could be defined as well. In 
dementia for example, the multi-state process might consider three states: healthy, demented and dead, 
and it would be necessary to consider the time spent in the demented state before death ( Commenges| 
et al. 20071. The joint multi-state model and its associated implemented function handle for semi- 


Markov. 


In summary, we introduce here a first joint model for longitudinal and multi-state clinical progression 
data. We showed that this model can easily be implemented under R and can be applied in practice 
through an example, the prostate cancer progression, which is one of many biomedical areas in which 
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such data are collected. This model that captures the complete information about the progression opens 
to much more precise knowledge of diseases and specific dynamic predictions. 


7 Appendix 

Appendix A: Simulation results 

In Section 4 in the article, the simulations were performed on 1500 subjects and showed that 9 Gauss- 
Hermite quadrature points should be used in this case. A complementary simulation study was based 
on the first 500 and 1000 individuals of each of the 500 replicates and used 9 Gauss-Hermite quadrature 
points. With 500 and 1000 subjects, we respectively observed 167 and 333 transitions 0^1, 103 and 
205 transitions 0 ^ 2, 96 and 191 transitions 1^2. The results are presented in Table|^ 

[Table 3 about here.] 

Overall, these results confirmed fhe good performances of fhe implemenfed funcfion. The coverage 
rafes were close fo 95% and fhe relafive bias were low, excepf for 702 ,x and ? 7 l 2 ,slope^ cerfainly due fo 
fhe required accuracy. 

Appendix B: Parametric versus non-parametric transition probabilities 

To assess fhe resulfs obfained in our application (Secfion 5 in fhe article), we compared fhe paramefric 
esfimafor of fhe fransifion probabilifies fo a non-paramefric estimator. The following resulfs are based 
on fhe book of ?. 


Preliminaries 


Consider fhe mulfi-sfafe process E = {E{t), t > 0) wifh values in fhe finite space S = [0,1,..., M} and 

where E{t) denotes fhe sfafe occupied by an individual af lime t. We assume fhaf E is a non-homogeneous 

Markov process, wifh leff Iruncafion and righf censoring. In fhe following, we will consider fhaf all fhe 

inlroduced mulfi-sfafe processes guaranlee fhe above properties. The inlensily of fransifion from sfafe 

Vv{E{t + dt) - k\E{t) - h) 


h € S to sfafe k e 5 al fime t is defined as Ahkit) - limdi 


t^o 


dt 


and we wrile 


/l(5, t) = {Ahk{s,t)} the (M -I- 1) X {M + 1) matrix of transition intensities. The matrix of cumulative 
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transition intensities is noted A(t), composed of non-diagonal elements Ahk{t) = Jg -Ihkiu) dn, V/i k, 
and diagonal elements Ahh(t) = - Yjki^h ^hk{t). Let us consider the transition probability Phki^, t) - 
Pr{E(t) - k\E{s) = h), with s <t, which is the probability that a subject in state h at time s occupies the 
state ^ at a later time t. We call P(5, t) - {Phki^, 0) the matrix of transition probabilities, which satisfies 
the Chapman-Kolmogorov equation: 

P(5, t) = P(5, m)P(m, t), with 0 < s < u < t. (8) 


Thus is deduced that P(5, t) is the unique solution of the Kolmogorov forward differential equations: 


Non-parametric estimator 


d 

dt 


P(5, s) ^ I, 

P(5, t) = P(5, t)A{t). 


(9) 


Let Nhk{t) be the number of direct observed transitions from state h to state k up to time t, and T/,(t) 
the number of individuals in state h just before time t. The non-parametric estimator of the cumulative 
intensities, called A*{t) has elements A*^^{t) estimated through the Nelson-Aalen estimator: 


Am(0 - 


Jo Yh{u) 


( 10 ) 


and = - Y.k*i, A’^^(f). The solution to the Kolmogorov equations (|^ permits to express the non- 

parametric estimate of the transition probabilities P*(i', t) though the product-integral: 


P*(5,0 = J((l + dA»). (11) 

where A*{u) is the non-parametric estimate of the cumulative transition intensities at time u, with 
dA^^(M) > -1 for all u, and I is the (M -i- 1) x (M + 1) identity matrix. This estimated matrix of 
transition probabilities is called Aalen-Johansen estimator. 

Let s < T* < ... < < t be the ordained times of observed direct transitions between s and t for 
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all individuals. From ([TT|), it can be deduced: 


/=! 


( 12 ) 


where AA*(T*) - A*{T*) - A*(r;_j) and AA^^iT*) > -1 for all T*. 
In our application, we applied: 


anoXt;) 
^ Fo(r;) 

0 


I + AA*{T*) - 


0 


0 

0 


anoi{t;) 

YoiT*) 

AAiXr;) 

Fi(r;) 

0 

0 

0 


ano2{t;) 

ano3{t;) 

ANoAiT;)\ 

Yo{t;) 

anuIt;) 

YoiT*) 

AAi3(r;) 

YoiT;) 

ANuiT;) 

Ti(r;) 

AfV2.(r;) 

YiiT*) 

aa23(f;) 

YliT;) 

AN2AiT;) 

Y2{t;) 

0 

YiiT;) 

^ AN3.iT;) 
Y3iT;) 

Y2iT;) 

AN34iT;) 

Y3iT;) 

0 

0 

1 


where Nh.{T*) - Y^u^hNhkiT;) and ANhkiT^) - NnkiT^) - Nhk{Tl^). 


Parametric estimator 

For each patient i e {1,..., A), the observed multi-state process is {£'^(0^ I ^ 0), where Ej(t) denotes the 
state occupied by subject i at time t and takes values in the finite space S = {0,1,..., M}. In our joint 
multi-state model, we were interested in the transitions intensities: 


= Jim 


Pr{Eiit + dt) = k\Eiit) = h\ bd 
dl ’ 


(13) 


which share the random effects bi with the longitudinal sub-model. Based on the Preliminaries para¬ 
graph, let A'{t\bi) = {A‘j^i^{t\bi)} be the parametric matrix of cumulative intensities for subject i, and 
P‘j^i^(s,t\bi) = {P‘f^f^{s,t\bi)} be the parametric matrix of transition probabilities. From the individual co¬ 
variates measured at baseline . and the parameters estimated by maximum likelihood 6, we computed 
for each subject i the individual predictions of the transition intensities and deduced the individ¬ 

ual predictions of the cumulative intensities Ajj^(t|6). To obtain the parametric estimator of the individual 
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transition probabilities P'(5, l\9), we could use the relations: 


P'hhi^, t\9) ^ exp - A^^(5|0)) , 

Pi,is,t\9) = J ^ k. 


(14) 


In practice, when the state space S is large, these integrals become too complex numerically. In the 
application we therefore preferred to calculate P'Ci', t\9) through the product-integral: 


P'(5, t|0) - J[ (l + dA'(M)), (15) 

(^,f] 

with > -1 for all u. The matrix of transition probabilities P(5, t) was then deduced by averaging 

over the N individual predictions: 

V{s,t\9) =-Yj^\s,t\h (16) 

1=1 

Covariance estimation of the non-parametric estimator 

The covariance matrix of the Aalen-Johansen estimator of transition probabilities can be estimated by 
the Greenwood-type estimator: 

COV(P*(5, t)) = P*(u, 0^ ® P*(5, M-)COV(dA*(M))P*(M, t)) ® P*(i', U-y, (17) 

where ® denotes the Kronecker product and ^ denotes the vector transpose. 

? simplified the computations in ([TT]) using the recursion formula: 


cov(P*(5,0) = {(I + AA*(t))^ ® I)cov(P*(5, ?-)){(! + AA*(0) ® I) + 
{I g) P*(5, t-))cov(AA*(0){I ® P*(5, t-)}. 


( 18 ) 
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where 


( {Yh{t) - ANhXmNh.it) 


eov(AA^^(0,AA;*,^,(0) ^ 


Yhitf 

{Yhit) - ANhimNhAt) 


, forh-k-h' = k'. 


Yhit? 

(Skk'Yhit) - ANhkit))ANhAt) 


, fork ^k^h' k'. 


Yhitf 


, for h = h' ,h k,h k', 


0, for h + h'. 


with 6kk' the Kroneeker delta. Note that cov(P*(5, t)) and cov(AA*(t)) are two {K + 1)^ x + 1)^ 
eovarianee matriees. 

These results may be used to eonstruet the 95% pointwise eonfidenee intervals for the Aalen- 
Johansen estimator: 


exp 


log 0 + 1.96 


V 




Appendix C: Example of R Code 


This example is based on the simulation study (see Seetion 4 in the artiele), with a non-homogeneous 
Markov multi-state model whieh ineluded three states and three transitions. The same eovariate, ealled 
X in the above eode, impaeted the longitudinal and survival proeesses, the longitudinal sub-model had 
a random intereept and a random slope, the log baseline intensities were approximated using B-splines, 
and the dependenee between the two proeesses was explained through the true eurrent level and the true 
eurrent slope of the biomarker. 

The JMstateModelO funetion and several detailed examples are available on https://github 
.com/LoicFerrer/JMstateModel/. 

# Load the packages and the function to estimate joint multi-state models: 
library(mstate) 

library(JM) 

sourceC"JMstateModel.R") 

# Import two databases which contain longitudinal and survival data: 
loadC'data. RData") 


############################### 
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#### Longitudinal sub-part #### 
############################### 

# Linear mixed model: 

ImeFit <- ImeCfixed = Y ~ (1 + times) •' X, 
data = data_long, 
random = ~ (1 + times) | id, 
method ^ "REML", 
control = listCopt = "optim")) 


############################## 

#### Multi-state sub-part #### 

############################## 

# Construct the 3*3 matrix of possible transitions: 
tmat <- matrixCNA, 3, 3) 

tmat[l, 2:3] <- 1:2 
tmat[2, 3] <- 3 

dimnames(tmat) <- list(from = c("State 0", "State 1", "State 2"), 

to = cC'State 0", "State 1", "State 2")) 

tmat 

# The transition ’Q -> 1’ is called ’1’,’® -> 2’ is called ’2’ and 

# ’1 -> 2’ is called ’3’. 

# Define the covariate(s) in the multi-state sub-part: 
covs <- "X" 

# The ’msprepO’ function divides the survival database in order to have 

# one line per transition at risk for each subject, with ’Tstart’ the 

# entry time in the current state, and ’Tstop’ the time of transition or 

# censorship; ’status’ denotes if the transition has been performed: 
data_mstate <- msprep(time = c(NA, "time_of_state_l", "time_of_state_2"), 

status = c(NA, "State 1", "State 2"), 

data = data_surv, 

trans = tmat, 

keep = covs, 

id = "id") 

# ’expand.covsC)’ permits to define the set of covariates which impacts 

# each transition: 

data_mstate <- expand.covs(data_mstate, covs, 

append = TRUE, longnames = FALSE) 

# Multi-state model with proportional hazards: 
coxFit <- coxph(Surv(Tstart, Tstop, status) 
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X.l + X.2 + X.3 + strata(trans), 
data = data_mstate, 
method = "breslow", 

X = TRUE, 
model = TRUE) 


#################################### 

#### Joint multi-state sub-part #### 

#################################### 

# Define the derivative of the fixed and random parts in the mixed model, 

# and indicate which covariates are kept, for the dependency on the slope 

# of the marker: 

dForm <- list(fixed = ~ 1 + X, 

indFixed = c(2, 4), 
random = ~ 1, 
indRandom =2) 

# Joint multi-state model with: 

# - current level and current slope as dependence function, 

# - adaptative cubic B-splines to approximate the log-baseline intensities, 

# - 9 Gauss-Hermite quadrature points in the pseudo-adaptative numerical 

# integration to approximate the integrals over random effects, 

# - 3 Gauss-Kronrod quadrature points to approximate the integral over time. 
jointFit <- 

JMstateModel(lmeObject = ImeFit, 
survObject = coxFit, 
timeVar = "times", 
parameterization = "both", 
method = "spline-PH-aGH", 

interFact = list(value = ~strata(trans) - 1, 
slope = ~strata(trans) - 1, 
data = data_mstate), 
derivForm = dForm, 

Mstate = TRUE, 
data.Mstate = data_mstate, 

ID.Mstate = "id", 

control = listCGHk = 9, Ing.in.kn =3)) 
summary(jointFit) 
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The funetion JMstateModelC), whieh is an extension of the standard jointModel() funetion to the 

multi-state framework, is available with several examples on https://github.com/LoicFerrGr/ 

JMstateModel/. 
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Table 1: Simulation results aeeording to 3, 9 and 15 pseudo-adaptative Gauss-Hermite quadrature points. For eaeh seenario, the statisties are 
(from left to right): mean, mean standard error, standard deviation, relative bias (in pereentage) and eoverage rate (in pereentage). 
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Figure 1: Simulated multi-state proeess. Arrows indieate the direetions of the possible transitions. Ahk{t) 
eharaeterizes the intensity of transition between states h and k at time t. The matrix Tgim has size (3,3) 
and is eomposed of elements Tsim,(ft+i)(t:+i), where Tsim,(/i+i)(<:+i) is the average number of observed 
direct transitions h ^ k over the 500 replicates. The diagonal elements Tsim,(/i+i)(ft+i) denote the average 
number of patients who were censored in state h. Note that the sum of elements of a row (/i -i-1) of T^m 
corresponds to the average number of patients who experienced the state h. 
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Figure 3: Multi-state representation of the elinieal progressions in prostate eaneer. Arrows indieate the 
directions of the possible transitions {N = 1474). Ahk{t) characterizes the intensity of transition between 
states h and k at time t. The matrix T has size (5,5) and is composed of elements T(/,+i)(,t+i)? where 
ft/i+iKAr+i) is the number of observed direct transitions h ^ k. The diagonal elements y(h+\)(h+\) denote 
the number of patients who were censored in state h. Note that the sum of elements of a row (/i -i- 1) of 
T corresponds to the number of patients who experienced the state h. 
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(a) Conditional standardized residuals versus the fitted val¬ 
ues. The black curve indicates the mean of conditional 
residuals according to the fitted values, using a locally 
weighted polynomial regression. 



Years since the end of EBRT 


(b) Observed and predicted values of the biomarker. 95% 
confidence intervals of the observed values are connected 
in grey. 
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(c) Predicted transition probabilities from the joint multi-state model and non-parametric probability transitions. 


Figure 4: Goodness-of-fit plots for the longitudinal process (a,b) and the multi-state process (c). 
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Table 2: Description of the two cohorts. 


Cohort 

RTOG 9406 

BCCA 

Pooled 

Study period 

1994-2013 

1994-2012 


Number of patients 

629 

845 

1474 

Number of PSA measures per patient 

13 (4, 23) 

9 (3, 15) 

10(3,21) 

iPSA* 

2.0 (1.0, 3.0) 

2.1 (0.6, 3.3) 

2.1 (0.8, 3.1) 

Clinical T-stage 

1 

355 (56.4%) 

184 (21.8%) 

539 (36.6%) 

2 

261 (41.5%) 

514 (60.8%) 

775 (52.6%) 

3^ 

13 (2.1%) 

147 (17.4%) 

160 (10.9%) 

Gleason score 

2-6 

424 (67.4%) 

605 (71.6%) 

1029 (69.8%) 

7 

167 (26.6%) 

189 (22.4%) 

356 (24.2%) 

8-10 

38 (6.0%) 

51 (6.0%) 

89 (6.0%) 

Mean time of first event^ 

9.8 (2.3, 15.9) 

7.7 (1.9, 14.1) 

8.2 (2.0, 15.0) 

Mean time of last contact* 

11.6 (2.9, 16.7) 

9.0 (3.4, 14.8) 

9.7 (3.1, 15.9) 


Continuous data: Median (5th and 95th percentiles). 

Categorical data: Amount (percentage). 

Times are in years since the end of EBRT. 

• Pre-therapy PSA value (ng/ml) in the log(. -i- 0.1) scale. 

^ Minimum between the time of first transition and the time of censoring. 
^ Minimum between the time of death and the time of censorship. 
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Table 3: Parameter estimates, standard errors and values in the joint multi-state model on the pooled 
data (N - 1474). 


Longitudinal Proeess Multi-state Proeess 



Value 

StdErr 

p-value 


Value 

StdErr 

p-value 

ySo 

-0.26 

0.06 

< 0.001 

'y02,iPSA 

0.34 

0.08 

< 0.001 

^O.iPSA 

0.80 

0.03 

< 0.001 

T04,iPSA 

0.27 

0.08 

< 0.001 

^O.cohort 

-0.01 

0.02 

0.525 

7(13,14,23,24,34),iPSA 

-0.24 

0.08 

0.002 

Pi 

0.71 

0.13 

< 0.001 

7(01,03,23),tstage2 

0.91 

0.18 

< 0.001 

^l.iPSA 

0.89 

0.06 

< 0.001 

7(01,03,23),tstage3-4 

0.73 

0.23 

0.002 

^l,tstage2 

0.38 

0.08 

< 0.001 

7(12,14,34),tstage2 

-0.04 

0.25 

0.865 

^l,tstage3-4 

0.48 

0.13 

< 0.001 

7(12,14,34),tstage3-4 

0.38 

0.30 

0.204 

^1,cohort 

-0.04 

0.04 

0.333 

7(03,23)gleason7 

0.97 

0.24 

< 0.001 

Pi 

-0.24 

0.04 

< 0.001 

7(03,23)gleason8- IS 

0.09 

0.43 

0.827 

^2,iPSA 

0.19 

0.02 

< 0.001 

7(01,14,24,34),cohort 

-0.42 

0.06 

< 0.001 

^2,tstage2 

0.14 

0.02 

< 0.001 

7(13,23),cohort 

0.89 

0.17 

< 0.001 

^2,tstage3-4 

0.26 

0.04 

< 0.001 

^(12,13) 

4.18 

0.38 

< 0.001 

^2,gleason7 

0.07 

0.02 

< 0.001 

^23 

3.03 

0.52 

< 0.001 

^2,gleason8-10 

0.22 

0.04 

< 0.001 

701,level 

0.37 

0.09 

< 0.001 

^2,cohort 

-0.06 

0.01 

< 0.001 

702,level 

0.51 

0.07 

< 0.001 

log(o-) 

-1.30 

0.01 


703,level 

0.45 

0.11 

< 0.001 





704,level 

-0.17 

0.05 

0.001 

Dn 

0.37 

0.02 


17l2,level 

-0.16 

0.10 

0.110 

Di2 

0.01 

0.01 


7l3,level 

-0.41 

0.20 

0.042 

T>13 

0.35 

0.03 


17l4,level 

0.10 

0.14 

0.487 

D22 

0.14 

0.01 


723,level 

-0.15 

0.09 

0.120 

D23 

0.25 

0.02 


724,level 

0.04 

0.05 

0.412 

D33 

1.68 

0.09 


734,level 

0.04 

0.08 

0.609 





701,slope 

2.54 

0.31 

< 0.001 





702,slope 

3.04 

0.25 

< 0.001 





703,slope 

2.43 

0.49 

< 0.001 





704,slope 

1.03 

0.32 

0.001 





7l2,slope 

2.01 

0.61 

0.001 





713,slope 

3.18 

0.80 

< 0.001 





7l4,slope 

-0.20 

1.27 

0.873 





1723,slope 

0.97 

0.67 

0.150 





724,slope 

0.29 

0.52 

0.583 





734,slope 

-0.79 

0.78 

0.313 


Dij denotes the /^'-element of the eovarianee matrix for the random effeets. y{hk,h'k'),x denotes the effeet 
of the eovariate X on the intensities of transitions h ^ k and h' k', i-e-y(hk,h'k'),x = Tm.x = yk'k'x- In 
the same idea, it is used ^(i 2 ,i 3 ) = ^i 2 - ^ 13 - 
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Table 4: Simulation results on the 500 and 1000 first individuals of the 500 replieates and using 9 Gauss- 
Hermite quadrature points. For eaeh seenario, the statisties depieted are, from left to right: mean, mean 
standard error, standard deviation, relative bias (in pereentage) and eoverage rate (in percentage). 


500 subjects 1000 subjects 



True 

value 

Mean 

StdErr 

StdDev 

Rel. 

bias 

Cov. 

rate 

Mean 

StdErr 

StdDev 

Rel. 

bias 

Cov. 

rate 

Longitudinal process 
ySo -0.793 

-0.797 

0.085 

0.086 

0.4 

94.8 

-0.798 

0.060 

0.061 

0.6 

94.4 

Po,x 

0.543 

0.544 

0.039 

0.040 

0.2 

94.8 

0.545 

0.028 

0.029 

0.3 

95.8 

Pi 

-0.096 

-0.096 

0.036 

0.036 

-0.2 

96.0 

-0.096 

0.026 

0.025 

0.1 

95.2 

Pl,X 

0.027 

0.027 

0.017 

0.017 

-0.1 

95.6 

0.027 

0.012 

0.012 

-0.1 

95.2 

log(o-) 

-0.737 

-0.737 

0.007 

0.019 

-0.1 

96.0 

-0.737 

0.005 

0.005 

-0.0 

93.0 

Multi-state 

7oi,x 

process 

0.281 

0.298 

0.136 

0.136 

5.9 

95.0 

0.296 

0.095 

0.094 

5.2 

95.2 

702,X 

0.023 

0.010 

0.155 

0.161 

-57.9 

93.0 

0.028 

0.108 

0.112 

21.6 

93.8 

712,X 

-0.169 

-0.162 

0.171 

0.176 

-4.3 

94.6 

-0.164 

0.117 

0.125 

-3.1 

94.8 

^01,level 

0.925 

0.927 

0.130 

0.144 

0.2 

93.0 

0.910 

0.090 

0.091 

-1.6 

94.6 

?702,level 

0.297 

0.298 

0.113 

0.109 

0.6 

96.0 

0.292 

0.079 

0.075 

-1.5 

95.6 

?712,level 

0.071 

0.071 

0.135 

0.141 

0.3 

94.4 

0.072 

0.092 

0.093 

0.4 

95.4 

^01,slope 

1.344 

1.514 

0.773 

0.827 

12.6 

93.6 

1.531 

0.541 

0.524 

13.9 

95.4 

^02,slope 

-1.096 

-1.071 

1.125 

1.119 

-2.3 

97.6 

-1.041 

0.788 

0.737 

-5.0 

95.6 

^12,slope 

0.009 

0.086 

1.425 

1.442 

853.3 

96.4 

0.036 

0.976 

1.007 

298.2 

94.6 

Random effects 

Dn 0.349 

0.346 

0.025 

0.025 

-0.8 

95.0 

0.347 

0.017 

0.017 

-0.6 

95.4 

D\2 

-0.041 

-0.042 

0.008 

0.008 

0.3 

95.0 

-0.042 

0.006 

0.006 

0.5 

95.8 

D 22 

0.062 

0.062 

0.004 

0.005 

-0.5 

93.8 

0.062 

0.003 

0.003 

-0.2 

92.6 



